% 衰变模拟器！（不是）
% 每时刻选取一定数量的未衰变原子，并使其发生衰变
% 虽然每个原子的衰变是独立的过程，但是统计上说，未衰变原子数越少，整体衰变速率越低，呈现一级反应特点（速率与原子个数成正比）。
% Gitee Repo

clc
clear

N = 50; %每行/列的原子个数，总共N*N个原子
TICK = 100; % 总运行步数
probability = 0.05; % 衰变概率（衰变速率）

[x y] = meshgrid(0:N-1);
A = ones(N,N);

figure()
hold on
axis equal
h = surf(x,y,A,'edgecolor','none');
caxis([0 1])
pause(1)

% 统计数据
% 列：时刻、衰变原子数、未衰变原子数
nums = zeros(TICK,3);
nums(:,1) = 1:TICK;

for tick = 1:TICK
  B = rand(N);
  C = A & B > (1-probability);
  A(C) = 0;

  % 未衰变原子=1、衰变原子=0
  set(h, 'ZData', A);
  set(h, 'CData', A);

  nums(tick,2) = sum(A(:)==0);
  nums(tick,3) = sum(A(:)==1);

  drawnow
  pause(0.1)

  if nums(tick,3) == 0
    nums(tick+1:end,2) = nums(tick,2);
    break;
  end
end

figure()
hold on
plot(nums(:,1),nums(:,3));
plot(nums(:,1),N*N*exp(-probability*nums(:,1)));
legend('模拟','理论（指数衰减）')
xlabel('时刻')
ylabel('未衰变原子个数')
ylim([0 N*N])

